import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.tsa.stattools import adfuller # ADF检验调用的命令
from pmdarima import auto_arima
import warnings
warnings.filterwarnings("ignore")
import matplotlib

# 设置全局字体
matplotlib.rcParams['font.sans-serif'] = ['SimHei']  # 使用黑体
matplotlib.rcParams['axes.unicode_minus'] = False  # 正确显示负号

# 设置随机种子以确保结果可重复
np.random.seed(42)

# === 1. 生成示例时间序列数据（无季节性） ===
dates = pd.date_range(start='2020-01-01', end='2024-12-01', freq='M')
n_periods = len(dates)

trend = np.linspace(100, 200, n_periods)
noise = np.random.normal(0, 5, n_periods)
data = trend + noise

ts_data = pd.DataFrame(data, index=dates, columns=['value'])
print("生成的时间序列数据预览：")
print(ts_data.head(10))

# 可视化原始数据
plt.figure(figsize=(12, 6))
plt.plot(ts_data, label='原始数据', color='blue')
plt.title('示例时间序列数据（无季节性）')
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()

# === 2. 检查时间序列的平稳性 ===
def adf_test(series, title=''):
    print(f'\n=== {title} ===')
    result = adfuller(series)
    print('ADF 统计量:', result[0])
    print('p 值:', result[1])
    print('临界值:', result[4])
    if result[1] <= 0.05:
        print("p 值 <= 0.05，序列是平稳的")
    else:
        print("p 值 > 0.05，序列是非平稳的，可能需要差分")

adf_test(ts_data['value'], '原始数据平稳性检验')
ts_diff = ts_data['value'].diff().dropna()
adf_test(ts_diff, '一阶差分后数据的平稳性检验')

plt.figure(figsize=(12, 6))
plt.plot(ts_diff, label='一阶差分数据', color='green')
plt.title('一阶差分后的时间序列')
plt.xlabel('日期')
plt.ylabel('差分值')
plt.legend()
plt.grid(True)
plt.show()

# === 3. 寻找最优 ARIMA 参数（基于 AIC） ===
print("\n=== 寻找最优 ARIMA 参数（基于 AIC） ===")

# 定义最优 AIC 和相应的最佳参数
best_aic = float('inf')  # 初始化为一个非常大的值
best_order = None
best_model_fit = None  # 保存最佳模型的拟合结果

# 遍历不同的 p 和 q 值
for p in range(0, 5):  # p 从 0 到 4
    for q in range(0, 5):  # q 从 0 到 4
        try:
            # 尝试构建 ARIMA 模型
            model = ARIMA(ts_data['value'], order=(p, 1, q))
            model_fit = model.fit()

            # 获取当前模型的 AIC 值
            current_aic = model_fit.aic
            print(f"ARIMA({p}, 1, {q}) 的 AIC: {current_aic:.2f}")

            # 如果当前 AIC 值更小，则更新最佳参数和模型
            if current_aic < best_aic:
                best_aic = current_aic
                best_order = (p, 1, q)
                best_model_fit = model_fit  # 保存最佳模型的拟合结果

        except Exception as e:
            print(f"ARIMA({p}, 1, {q}) 模型拟合失败：{e}")
            continue

print(f"\n最佳 ARIMA 参数（基于 AIC）：{best_order}，AIC: {best_aic:.2f}")

# === 4. 使用最佳参数拟合 ARIMA 模型 ===
# 直接使用保存的最佳模型拟合结果
model = ARIMA(ts_data['value'], order=best_order)
model_fit = best_model_fit

print("\n=== 最佳 ARIMA 模型摘要 ===")
print(model_fit.summary())

# === 5. 预测未来 12 个月（含置信区间） ===
forecast_result = model_fit.get_forecast(steps=12)
forecast_mean = forecast_result.predicted_mean
forecast_ci = forecast_result.conf_int()

# 直接使用 forecast_mean 的索引
forecast_series = pd.Series(forecast_mean.values, index=forecast_mean.index)

plt.figure(figsize=(12, 6))
plt.plot(ts_data, label='历史数据', color='blue')
plt.plot(forecast_series, label='预测值', color='red', linestyle='--')
plt.fill_between(
    forecast_ci.index,
    forecast_ci.iloc[:, 0],
    forecast_ci.iloc[:, 1],
    color='red', alpha=0.1, label='95% 置信区间'
)
plt.title('最佳 ARIMA 模型预测（2025年，含置信区间）')
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()

print("\n=== 2025 年预测值 ===")
forecast_df = pd.DataFrame(forecast_series, columns=['预测值'])
print(forecast_df)

# === 6. 模型评估（使用训练/测试集分割，基于 MAPE 和其他指标） ===
train = ts_data['value'][:-12]
test = ts_data['value'][-12:]

model_train = ARIMA(train, order=best_order)
model_train_fit = model_train.fit()

forecast_test_result = model_train_fit.get_forecast(steps=12)
forecast_test_mean = forecast_test_result.predicted_mean
forecast_test_ci = forecast_test_result.conf_int()
forecast_test_series = pd.Series(forecast_test_mean.values, index=forecast_test_mean.index)

# 定义 calculate_mape 函数
def calculate_mape(actual, predicted):
    actual, predicted = np.array(actual), np.array(predicted)
    mask = actual != 0
    if np.sum(mask) == 0:
        return np.inf
    return np.mean(np.abs((actual[mask] - predicted[mask]) / actual[mask])) * 100

# 计算 MAPE、MAE 和 RMSE
from sklearn.metrics import mean_absolute_error, mean_squared_error

mape = calculate_mape(test, forecast_test_series)
mae = mean_absolute_error(test, forecast_test_series)
rmse = np.sqrt(mean_squared_error(test, forecast_test_series))

print(f"\n=== 模型评估 ===")
print(f"测试集 MAPE（平均绝对百分比误差）: {mape:.2f}%")
print(f"测试集 MAE（平均绝对误差）: {mae:.2f}")
print(f"测试集 RMSE（均方根误差）: {rmse:.2f}")

# 可视化测试集预测结果（含置信区间）
plt.figure(figsize=(12, 6))
plt.plot(train, label='训练数据', color='blue')
plt.plot(test, label='测试数据', color='green')
plt.plot(forecast_test_series, label='测试集预测值', color='red', linestyle='--')
plt.fill_between(
    forecast_test_ci.index,
    forecast_test_ci.iloc[:, 0],
    forecast_test_ci.iloc[:, 1],
    color='red', alpha=0.1, label='95% 置信区间'
)
plt.title('ARIMA 模型测试集预测效果（含置信区间）')
plt.xlabel('日期')
plt.ylabel('值')
plt.legend()
plt.grid(True)
plt.show()

# === 7. 残差分析 ===
residuals = model_train_fit.resid

plt.figure(figsize=(12, 6))
plt.plot(residuals, label='残差')
plt.title('ARIMA 模型残差')
plt.xlabel('日期')
plt.ylabel('残差值')
plt.legend()
plt.grid(True)
plt.show()

from statsmodels.graphics.tsaplots import plot_acf
plot_acf(residuals, lags=20)
plt.title('残差自相关图')
plt.show()